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1.  Objectives: 


The  main  objective  of  this  research  was  to  develop  and  apply  high  order  accuracy  methods 
(Spectral  and  Finite  Differences)  to  the  numerical  simulation  of  flows  with  discontinuities,  in 
complex  geometries.  The  prime  target  of  this  effort  is  the  simulations  of  supersonic  reactive 
flows. 

The  work  entailed  research  in  different  areas,  the  results  of  these  converge  in'  the  codes 
developed.  Some  of  these  areas  are: 

•  The  application  of  spectral  methods  to  shock  wave  calculations. 

•  Development  of  high  order  accuracy  finite  differences  ENO  and  YVENO  schemes. 

•  The  development  of  the  Discontinuous  Galerkin  methods. 

•  The  resolution  of  the  Gibbs  Phenomenon. 

•  The  developments  of  high  order  far  fieid  boundary  conditions  and  in  particular  ab¬ 
sorbing  layers  for  aero-acoustic  problems. 

•  General  research  into  the  theory  of  approximations  of  PDE’s. 

The  above  research  yielded  results  that  were  applied  in  different  problems  relevant  to 
US  Air-force.  As  an  example  we  show  how  to  get  rid  of  the  oscillations  in  picture  splicing 
(see  later). 


2.  Summary  of  Research: 


•  Spectral  Simulations  of  Supersonic  Reactive  Flows 

The  culmination  of  the  research  effort  under  this  grant  is  the  construction  of  a  multidi¬ 
mensional  spectral  code  for  the  simulations  of  complicated  interactions  of  shock  waves 
and  reactive  flows.  A  three  dimensional  supersonic  reactive  Navier-Stokes  Solver  using 
Chebyshev  collocation  methods  has  been  used  in  the  study  of  mixing  in  a  Scramjet 
engine  and  addressed  the  issue  of  mixing  enhancement  by  shock  interactions.  The 
code  runs  on  the  IBM-SP2  parallel  computer  and  its  accuracy  and  stability  have  been 
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verified.  It  had  been  shown  that  the  code  provided  superior  information  to  low  order 
finite  difference  schemes. 

A  new  mechanism  that  might  be  responsible  for  breaking  up  the  fuel  integrity  (enhanc¬ 
ing  the  mixing)  was  identified.  The  heavier  fluids  ( H2O  and  02)  are  accelerated  by  the 
shock  penetrating  into  the  lighter  fluid  (H2)  and  tends  to  form  fingers  (related  to  the 
Richtmyer-Meshkov  instability).  As  the  flame  jet  reaches  the  other  side  of  the  fuel-air 
interface,  a  roll-up  vortex  (related  to  the  Kevin-Helmholtz  instability)  is  formed  at  the 
tip  of  the  flame  jet.  The  flame  jet  cuts  through  the  hydrogen  fuel  and  makes  contact 
with  the  air  on  the  other  side  of  the  fuel-air  interface.  The  vortex  at  the  tip  of  the 
flame  jet  lifts  and  breaks  up  the  fuel-air  interface.  The  motion  of  the  vortex  creates 
another  flame  jet  allowing  fresh  air  to  penetrate  inside  and  to  interact  with  the  fuel. 
This  process  repeats  itself  as  the  new  flame  jet  reaches  the  other  side  of  the  fuel-air 
interface. 

Low  order  schemes  tend  to  suppress  the  formation  of  the  vortex  due  to  its  inherently 
large  numerical  dissipation.  Therefore,  the  low  order  scheme  predicts  the  large  scale 
vortical  roll-up  but  not  the  break-up  of  the  fuel. 

The  results  of  this  research  were  reported  in  [9]  and  [10]. 

The  filtering  techniques  developed  for  this  problem  has  a  wide  range  of  applications  . 
In  fact  in  [42]  a  class  of  filters  based  upon  the  numerical  solution  of  high-order  elliptic 
problems  in  Rd  which  allow  for  independent  determination  of  order  and  cut-off  wave 
number  and  which  default  to  classical  Fourier-based  filters  in  homogeneous  domains 
were  applied.  These  filters  are  not  restricted  to  applications  in  tensor-product  based 
geometries  as  is  generally  the  case  for  Fourier-based  filters.  The  discrete  representation 
of  the  filtered  output  is  constructed  from  a  Krylov  space  generated  in  solving  a  well- 
conditioned  system  arising  from  a  low-order  PDE. 

•  Spectral  Methods  for  Complex  Geometries: 

A  major  difficulty  in  the  application  of  high  order  methods  to  realistic  problems  is  the 
issue  of  applying  high-order  formulations  to  complex  geometries.  Often,  generating 
a  reasonable  grid  around  a  complex  configuration  is  the  most  difficult  aspect  of  the 
solution  procedure. 

A  multi-domain  approach  based  on  breaking  the  geometry  into  piecewise  smooth  “sub- 
domains”,  using  quadrilaterals  (hexahedron)  in  two  (three)  dimensions  was  developed. 
Each  sub-domain  is  then  discretized  with  a  stable  tensor  product  formulation,  and  the 
resulting  sub-domains  are  patched  together. 
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It  turns  out  that  the  most  important  ingredient  of  this  method  is  how  to  impose  inter¬ 
face  boundary  conditions.  This  is  especially  important  for  high  order  schemes  as  an 
improper  imposition  of  those  boundary  conditions  can  lead  to  reduction  of  accuracy  as 
well  as  time  instabilities.  The  issue  is  very  subtle  as  the  scheme  may  still  be  classically 
stable,  but  display  non-physical  growth  in  time.  We  have  developed  a  methodology 
for  a  time  stable  and  accurate  imposition  of  interface  boundary  conditions.  It  is  valid 
for  high-order  finite-difference  (FD)  discretizations  and  certain  spectral  formulations. 
Our  method  is  based  on  a  penalty  formulation  where  the  penalty  parameters  are  de¬ 
termined  by  stability  considerations  or  other  properties  of  the  numerical  scheme.  The 
SAT  procedure  assures  time  stability  for  systems  of  equations  that  have  a  bounded 
energy  norm.  This  is  not  true  in  general  for  other  high-order  FD  methods.  Indeed,  non¬ 
penalty  approaches  often  lead  to  non-physical  growth  in  time  for  systems  of  equations, 
even  though  the  discretization  operator  is  stable  for  the  scalar  case.  The  situation  is 
even  more  serious  for  calculations  of  flows  with  shock  waves.  Currently,  there  are  no 
good  methods  that  can  pass  shock  waves  from  one  sub-domain  to  the  other  within 
high  order  schemes.  A  progress  had  been  made  for  this  problem,  and  at  least  for  the 
range  of  strength  of  shocks  relevant  to  reactive  flows,  a  procedure  was  developed.  For 
finite  difference  schemes  we  had  outlined  a  stable  and  conservative  interface  treatment 
of  arbitrary  spatial  accuracy. 

In  [3]  a  method  to  construct  Spectral  methods  for  arbitrary  grids  was  introduced, 
this  extends  the  validity  of  the  Spectral  methods  to  domains  of  arbitrary  shape.  In 
[32]  an  optimal  set  of  points  for  interpolation  in  triangles  were  found  and  using  the 
methodology  in  [3],  an  efficient  spectral  method  for  triangles  was  developed  in  [35]. 

This  advance,  together  with  the  multi-domain  methodology  described  above  extends 
spectral  methods  to  complicated  geometries. 

In  [24]  the  optimal  strategy  of  subdivision  of  domains  for  spectral  calculations  of  wave 
phenomena  was  discussed.  A  formula  connecting  the  optimal  number  of  sub-domains 
with  the  complexity  of  the  problem  (number  of  waves)  and  the  required  accuracy 
was  given.  In  [20]  this  was  extended  to  parallel  computers,  taking  the  number  of 
processors  and  communication  time  into  account.  It  was  shown  that,  for  present  day 
multicomputers,  the  impact  of  communication  overhead  does  not  significantly  shift 
the  number  of  domains  and  the  points  at  each  domain  from  the  optimal  uni-processor 
values,  and  that  the  effects  of  granularity  are  more  important.  A  different  approach 
to  multi-domain  methods  is  considered  in  [36]  where  a  wavelet  optimized  adaptive 
multi-domain  method  had  been  presented. 
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In  [28,  29,  33,  34]  the  penalty  method  for  the  multi-domain  interface  boundary  con¬ 
ditions  had  been  presented  for  the  Navier  Stokes  equations  as  well  as  for  general 
hyperbolic  equations  as  the  Maxwell’s  equations  for  electro-magnetics. 

ENO  and  Weighted  ENO  Schemes  and  Related  Topics 

We  have  performed  research  on  high  order  finite  difference  and  finite  volume  ENO 
(essentially  non-oscillatory)  and  WENO  (weighted  essentially  non-oscillatory)  schemes. 
These  are  schemes  suitable  for  the  computation  of  solutions  containing  both  shocks 
and  other  discontinuities  and  detailed  smooth  structures. 

ENO  idea  is  an  adaptive  interpolation  procedure  which  tries  to  automatically  choose  a 
locally  smoothest  region  to  perform  a  high  order  interpolation,  hence  avoiding 'crossing 
a  discontinuity  whenever  possible.  WENO  is  a  modification  and  improvement  of  ENO 
schemes.  Instead  of  using  only  one  of  the  many  candidate  stencils  based  on  local 
smoothness  as  in  ENO,  WENO  uses  a  linear  combination  of  the  contribution  from  all 
candidate  stencils,  each  with  suitable  nonlinear  weight. 

In  [46],  Shu  and  Zeng  have  applied  ENO  method  to  the  viscoelastic  model  with  fading 
memory.  The  memory  term  is  treated  by  introducing  new  variables  and  rewrite  the 
system  by  adding  more  differential  equations  but  without  explicit  memory  terms.  The 
appearance  of  the  memory  terms  regularizes  the  solution  somewhat,  and  in  many  cases 
it  is  still  a  theoretically  open  question  whether  shocks  will  develop  from  smooth  initial 
data.  We  have  performed  theoretical  analysis  about  the  linearized  system  for  large 
time,  and  have  applied  ENO  scheme  to  study  the  nonlinear  system  for  both  local  time 
and  large  time.  The  high  order  accuracy  and  sharp,  non-oscillatory  shock  transition 
allow  us  to  obtain  fine  resolution  for  tens  of  thousands  of  time  steps,  and  to  study  the 
shock  interactions  after  the  formation  of  shocks. 

Application  of  ENO  scheme  to  the  study  of  shock  longitudinal  vortex  interaction 
problem  is  carried  out  by  Erlebacher,  Hussaini  and  Shu  in  [13].  We  have  studied 
the  shock/longitudinal  vortex  interaction  problem  in  axisymmetric  geometry.  Linear 
analysis,  shock  fitting  code,  and  shock  capturing  ENO  are  used  in  different  parameter 
range,  to  study  various  cases  of  nearly  linear  regime,  weakly  nonlinear  regime,  and 
strong  nonlinear  regime.  Vortex  breakdown  as  a  function  of  Mach  number  ranging 
from  1.3  to  10  is  studied,  extending  the  range  of  existing  results.  For  vortex  strengths 
above  a  critical  value,  a  triple  point  forms  on  the  shock,  leading  to  a  Mach  disk.  This 
leads  to  a  strong  recirculating  region  downstream  of  the  shock.  It  is  found  out  that  a 
secondary  shock  forms,  to  provide  the  necessary  deceleration  so  that  the  fluid  velocity 


can  adjust  to  downstream  conditions  at  the  shock.  Also  on  ENO  schemes,  Harabetian, 
Osher  and  Shu  have  investigated  a  novel  Eulerian  approach  for  simulating  vortex  mo¬ 
tion  using  a  level  set  regularization  procedure  [27].  Our  approach  uses  a  decomposition 
of  the  vorticity  of  the  form  £  =  P{p)v-,  in  which  both  p  (the  level  set  function)  and 
T)  (the  vorticity  strength  vector)  are  smooth.  We  derive  coupled  equation.-,  for  p  and 
rj  which  give  a  regularization  of  the  problem.  The  regularization  is  topological  and 
is  automatically  accomplished  through  the  use  of  numerical  schemes  whose  viscosity 
shrinks  to  zero  with  grid  size.  There  is  no  need  for  explicit  filtering,  even  when  sin¬ 
gularities  appear  in  the  front.  The  method  also  has  the  advantage  of  automatically 
allowing  topological  changes  such  as  merging  of  surfaces.  Numerical  examples  includ¬ 
ing  two  and  three  dimensional  vortex  sheets,  two  dimensional  vortex  dipole  sheets  and 
point  vortices,  are  given.  To  our  knowledge,  this  is  the  first  three  dimensional  vortex 
sheet  calculation  in  which  the  sheet  evolution  feeds  back  to  the  calculation  of  the  fluid 
velocity. 

In  [40],  Jiang  and  Shu  have  investigated  WENO  (weighted  ENO)  schemes,  extending 
the  ideas  of  Liu,  Osher  and  Chan.  In  [40],  the  weights  are  chosen  so  that  in  smooth 
regions,  including  at  smooth  local  extrema,  they  are  close  to  an  optimal  linear  weight 
which  gives  the  highest  possible  order  of  accuracy  of  an  upwind-biased  linearly  stable 
scheme.  Near  shocks,  however,  those  stencils  which  contain  the  shock  are  assigned 
essentially  zero  weights.  Thus  WENO  resembles  a  linear  high  order  upwind  biased 
scheme  in  smooth  regions,  and  resembles  ENO  near  shocks,  with  a  smooth  numerical 
flux  function.  One  important  advantage  of  WENO,  due  to  its  smoothness  of  fluxes, 
is  that  convergence  for  smooth  solutions  can  be  proven.  Also,  convergence  towards 
steady  states  is  easier  than  ENO. 

Also  about  high  order  weighted  essentially  non-oscillatory  (WENO)  schemes,  jointly 
with  Philippe  Montarnal,  we  have  used  a  recently  developed  energy  relaxation  theory 
by  Coquel  and  Perthame  and  high  order  weighted  essentially  non-oscillatory  (WENO) 
schemes  to  simulate  the  Euler  equations  of  real  gas  [43].  The  main  idea  is  an  energy 
decomposition  under  the  form  e  =  e\  -f  e2 ,  where  £x  is  associated  with  a  simpler 
pressure  law  (7-law  in  our  case)  and  the  nonlinear  deviation  £2  is  convected  with  the 
flow. 

A  relaxation  process  is  performed  for  each  time  step  to  ensure  that  the  original  pressure 
law  is  satisfied.  The  necessary  characteristic  decomposition  for  the  high  order  WENO 
schemes  is  performed  on  the  characteristic  fields  based  on  the  £1  7-law.  The  algorithm 
only  calls  for  the  original  pressure  law  once  per  grid  point  per  time  step,  without  the 
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need  to  compute  its  derivatives  or  any  Riemann  solvers.  Both  one  and  two  dimensional 
numerical  examples  are  shown  to  illustrate  the  effectiveness  of  this  approach. 

About  high  order  weighted  essentially  non-oscillatory  (WENO)  finite  volume  schemes 
on  general  triangulations,  Hu  and  Shu  [37],  [38]  constructed  third  and  fourth  order 
WENO  schemes  on  two  dimensional  unstructured  meshes  (triangles)  in  the  finite  vol¬ 
ume  formulation.  The  third  order  schemes  are  based  on  a  combination  of  linear  polyno¬ 
mials  with  nonlinear  weights,  and  the  fourth  order  schemes  are  based  on  combination 
of  quadratic  polynomials  with  nonlinear  weights.  We  have  addressed  several  difficult 
issues  associated  with  high  order  WENO  schemes  on  unstructured  mesh,  including  the 
choice  of  linear  and  nonlinear  weights,  grouping  techniques  to  avoid  negative  weights, 
etc.  Numerical  examples  are  shown  to  demonstrate  the  accuracies  and  robustness  of 
the  methods  for  shock  calculations. 

As  a  related  topic,  In  [26],  S.  Gottlieb  and  Shu  further  explored  a  class  of  high  or¬ 
der  TVD  (total  variation  diminishing)  Runge-Kutta  time  discretization  suitable  for 
solving  hyperbolic  conservation  laws  with  stable  spatial  discretizations.  We  illustrate 
with  numerical  examples  that  non-TVD  but  linearly  stable  Runge-Kutta  time  dis¬ 
cretization  can  generate  oscillations  even  for  TVD  (total  variation  diminishing)  spatial 
discretization,  verifying  the  claim  that  TVD  Runge-Kutta  methods  are  important  for 
such  applications.  We  then  explore  the  issue  of  optimal  TVD  Runge-Kutta  methods 
for  second,  third  and  fourth  order,  and  for  low  storage  Runge-Kutta  methods.  On 
another  related  topic,  Perthame  and  Shu  [44]  have  investigated  the  issue  of  positivity 
preserving  (for  density  and  pressure)  high  order  methods  for  compressible  Euler  equa¬ 
tions  of  gas  dynamics  on  arbitrary  triangulation.  A  general  framework  for  positivity 
is  established  and  examples  within  this  framework  are  given. 

•  Discontinuous  Galerkin  Methods 

A  quite  successful  technique  for  hyperbolic  conservation  laws  is  the  discontinuous 
Galerkin  finite  element  method.  In  this  method  the  partial  differential  equation  is 
multiplied  by  a  test  function,  integrated  over  a  cell,  and  formally  integrated  by  parts 
to  obtain  a  weak  formulation.  A  solution  is  sought  among  discontinuous  (across  cell 
interface)  piecewise  polynomials  of  r-th  degree  for  a  (r  +  l)-th  order  method.  Be¬ 
cause  of  the  discontinuity  at  cell  interface,  this  method  can  accommodate  successful 
finite  difference  methodology  (approximate  Riemann  solvers  and  limiters)  into  a  finite 
element  framework.  Theoretical  results  similar  to  finite  difference  methods,  such  as 
total  variation  stability  for  ID  and  maximum  norm  stability  for  2D  and  3D,  can  be 
proved  for  this  class  of  discontinuous  Galerkin  methods  of  arbitrary  order  of  accuracy 
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and  for  (almost)  arbitrary  triangulations.  An  essential  difference  between  this  class 
of  finite  element  method  and  the  finite  volume  method  (which  can  also  be  defined 
on  an  arbitrary  triangulation)  is  that  the  latter  has  only  one  independent  degree  of 
freedom  (the  cell  average)  over  each  cell,  while  the  former  has  many  (for  example,  it 
has  three  degrees  of  freedom  for  the  piecewise  linear  case  in  2D).  This  fact  renders 
the  scheme  more  local  (no  wide  stencil  reconstruction  is  needed  to  compute  the  flux 
at  cell  interface),  hence  more  suitable  for  parallel  computing,  and  provides  a  differ¬ 
ent  setting  for  theoretical  justification  of  stability  and  convergence  of  the  algorithm. 
In  practice,  finite  element  methods  can  handle  complicated  geometry  and  boundary 
conditions  more  easily. 

In  [6],  Cockburn  and  Shu  have  studied  the  Local  Discontinuous  Galerkin  methods 
for  nonlinear,  time-dependent  convection-diffusion  systems.  These  methods  are  an 
extension  of  the  Runge-Kutta  Discontinuous  Galerkin  methods  for  purely  hyperbolic 
systems  to  convection-diffusion  systems  and  share  with  those  methods  their  high  par- 
allelizability,  their  high-order  formal  accuracy,  and  their  easy  handling  of  complicated 
geometries,  for  convection  dominated  problems.  It  is  proven  that  for  scalar  equa¬ 
tions,  the  Local  Discontinuous  Galerkin  methods  are  L2-stable  in  the  nonlinear  case. 
Moreover,  in  the  linear  case,  it  is  shown  that  if  polynomials  of  degree  k  are  used, 
the  methods  are  &-th  order  accurate  for  general  triangulations;  although  this  order  of 
convergence  is  suboptimal,  it  is  sharp  for  the  LDG  methods.  Preliminary  numerical 
examples  displaying  the  performance  of  the  method  are  shown. 

In  [5],  Cockburn  and  Shu  have  extended  the  Runge-Kutta  discontinuous  Galerkin 
method  to  multidimensional  nonlinear  systems  of  conservation  laws.  The  algorithms 
are  described  and  discussed,  including  algorithm  formulation  and  practical  implemen¬ 
tation  issues  such  as  the  numerical  fluxes,  quadrature  rules,  degrees  of  freedom,  and 
the  slope  limiters,  both  in  the  triangular  and  the  rectangular  element  cases.  Numeri¬ 
cal  experiments  for  two  dimensional  Euler  equations  of  compressible  gas  dynamics  are 
presented  that  show  the  effect  of  the  (formal)  order  of  accuracy  and  the  use  of  triangles 
or  rectangles,  on  the  quality  of  the  approximation. 

In  [2],  Atkins  and  Shu  have  discussed  a  discontinuous  Galerkin  formulation  that  avoids 
the  use  of  discrete  quadrature  formulas.  The  application  is  carried  out  for  one  and  two 
dimensional  linear  and  nonlinear  test  problems.  This  approach  requires  less  compu¬ 
tational  time  and  storage  than  conventional  implementations  but  preserves  the  com¬ 
pactness  and  robustness  inherent  in  the  discontinuous  Galerkin  method. 

Hu  and  Shu  have  presented  a  discontinuous  Galerkin  finite  element  method  for  solving 
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the  nonlinear  Hamilton-Jacobi  equations  in  [39].  This  method  is  based  on  the  Runge- 
Kutta  discontinuous  Galerkin  finite  element  method  for  solving  conservation  laws.  The 
method  has  the  flexibility  of  treating  complicated  geometry  by  using  arbitrary  trian¬ 
gulation,  can  achieve  high  order  accuracy  with  a  local,  compact  stencil,  and  are  suited 
for  efficient  parallel  implementation.  One  and  two  dimensional  numerical  examples 
are  given  to  illustrate  the  capability  of  the  method.  In  [41],  Lepsky,  Hu  and  Shu  have 
further  investigated  this  method  from  theoretical  and  computational  points  of  view. 
Theoretical  results  on  accuracy  and  stability  properties  of  the  method  are  proven  for 
certain  cases  and  related  numerical  examples  are  presented.  It  should  be  noted  that 
for  spectral  methods,  the  penalty  imposition  of  boundary  conditions  is  identical  with 
the  DG  method. 

•  Psuedopack  -  Numerical  Library  for  Spectral  Differentiations: 

A  software  library  using  the  latest  and  best  algorithms  for  computing  Chebyshev,  Leg¬ 
endre  and  Fourier  derivative  for  multiple  data  set  with  optimal  accuracy  and  efficiency 
was  written.  This  is  important  since  spectral  methods  based  on  orthogonal  polynomial 
are  very  sensitive  to  roundoff  error.  Special  numerical  techniques  and  algorithms  were 
employed  to  increase  the  efficiency  and  accuracy  of  the  underlining  methods. 

The  package  has  the  following  features: 

1.  Fourier,  Chebyshev  and  Legendre  methods  on  the  Gauss-Lobatto  points  are 
supported. 

2.  Matrix-Matrix  Multiply  Algorithm,  Even-Odd  Decomposition  Algorithm  and 
Fast  Fourier/ Cosine  Transform  Algorithm  are  supported  for  computing  the 
derivative  of  a  function. 

3.  Compiled  on  IBM  RS/6000,  CRAY,  SGI,  SUN  and  Generic  UNIX  machine. 

4.  Native  fast  assembly  library  call,  when  available,  is  used  for  the  library’s 
computational  kernel. 

5.  Special  fast  algorithms  are  provided  for  cases  when  the  function  has  either 
even  or  odd  symmetry. 

6.  Mapping  was  used  to  reduce  the  roundoff  error  for  the  Chebyshev  and  Leg¬ 
endre  differentiation. 

7.  Extensive  built-in/User  definable  coordinate  transformation  routines. 

8.  Built-in  filtering  for  smoothing  of  the  function  and  its  derivative. 

9.  Unified  subroutine  call  interface  allows  modification  of  any  parameters  with¬ 
out  any  change  to  be  made  to  the  subroutine  call  statement. 
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10.  Simple  user  callable  subroutines  return  the  derivatives  of  a  multiple  data  set. 

Since  the  user  is  shielded  from  any  coding  errors  of  the  main  derivative  routines, 
reliability  of  the  solutions  is  enhanced.  It  speeds  up  code  development,  increases 
productivity  and  enhances  re-usability. 

The  package  is  available  at  www.cfm.brown.edu/people/wsdon/home.html.  An  early 
description  of  the  package  can  be  found  in  [12]. 

Some  of  the  ideas  that  were  incorporated  into  the  design  of  this  software  may  be  found 
in  [11]  where  accuracy  enhancement  for  higher  derivatives  using  Chebyshev  collocation 
and  a  mapping  technique  is  discussed.  Also  in  [7]  the  accurate  computations  of  high 
order  derivative  by  spectral  methods  (known  to  have  troubles  with  roundoff  errors) 
were  discussed. 

•  Acoustics 

Significant  progress  has  been  made  in  the  case  of  plane  acoustics  embedded  in  uniform 
flows.  First  we  tried  utilizing  the  well-behaved  PML  methods  from  ambient  acoustics  in 
combination  with  a  layer  that  slowly  decelerates  the  waves  to  a  zero  Mach  number  prior 
to  entering  the  PML  layer,  justifying  the  use  of  the  ambient  PML.  This  is  reported  in 
[30]  While  computations  confirm  the  efficiency  and  simplicity  of  the  proposed  scheme  it 
cannot  be  claimed  to  be  a  true  PML  method,  as  the  PML  property  is  obtained  only  for 
very  wide  layers.  We  then  applied  the  methodology  based  on  mathematical  method, 
to  develop  the  first  strongly  well-posed  PML  method  for  the  problem  of  advective 
acoustics.  The  additional  degrees  of  freedom,  required  to  ensure  the  PML  property, 
are  introduced  through  a  number  of  additional  ordinary  differential  equations,  and 
a  single  partial  differential  equation.  The  additional  equations  are  defined  so  as  to 
ensure  that  the  total  set  of  equations  support  decaying  wave  solutions,  independent  of 
frequency  and  angle  of  incidence  of  the  incoming  wave.  This  research  appears  in  [1]. 

It  is  clear  that  the  development  of  perfectly  matched  absorbing  layers  for  the  equa¬ 
tions  of  acoustics  is  in  its  infancy.  While  the  recent  encouraging  results  suggest  the 
possibility  of  developing  such  layers  for  a  variety  of  flow  conditions,  many  important 
questions  remain  open. 

The  work  drew  the  attention  of  researchers  from  Pratt  and  Whitney  who  are  working 
on  the  problem  of  accurate  modeling  of  turbine  flutter,  where  the  variable  mean  flow 
is  obtained  from  a  direct  solution  of  the  Euler  equations  and  the  noise  propagation 
problem  is  traced  in  this  mean  field  using  the  linearized  equations.  A  joint  work  is 
being  planned. 
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•  The  Resolution  of  the  Gibbs  Phenomenon 

The  nonuniform  convergence  of  the  Fourier  series  for  discontinuous  functions,  and  in 
particular  the  oscillatory  behavior  of  the  finite  sum,  was  already  analyzed  by  Wilbra- 
ham  in  1848.  This  was  later  named  The  Gibbs  Phenomenon. 

In  [22]  a  review  of  the  Gibbs  phenomenon  from  a  different  perspective  was  given. 
In  this  view,  the  Gibbs  phenomenon  deals  with  the  issue  of  recovering  point  values 
of  a  function  from  its  expansion  coefficients.  Alternatively  it  can  be  viewed  as  the 
possibility  of  the  recovery  of  local  information  from  global  information.  The  main 
theme  here  is  not  the  structure  of  the  Gibbs  oscillations  but  the  understanding  and 
resolution  of  the  phenomenon  in  a  general  setting. 

The  purpose  of  this  article  was  to  review  the  Gibbs  phenomenon  and  to  show  that 
the  knowledge  of  the  expansion  coefficients  is  sufficient  for  obtaining  the  point  values 
of  a  piecewise  smooth  function,  with  the  same  order  of  accuracy  as  in  the  smooth 
case.  This  is  done  by  using  the  finite  expansion  series  to  construct  a  different,  rapidly 
convergent,  approximation. 

In  [23]  a  general  method  of  extracting  high  quality  approximations  from  a  slowly  con¬ 
verging  ones  have  been  introduced.  Conditions  to  determine  which  low  order  approxi¬ 
mation  contains  enough  information  such  that  a  better  approximation  can  be  derived 
from  it  were  given.  Previous  results  on  the  resolution  of  the  Gibbs  phenomenon  were 
shown  to  be  special  cases  of  this  general  theory.  This  work  generalizes  and  extends 
the  previous  work  on  the  Gibbs  phenomenon,  and  the  supersonic  reactive  codes  use 
the  theory  to  overcome  the  Gibbs  oscillation  from  the  shock. 

Another  application,  motivated  by  a  request  from  R.  Albanese  (MPD  Brooks  AFB) 
of  the  problem  of  ’’gluing”  of  spliced  pictures  was  considered.  A  typical  situation  here 
is  that  a  two  dimensional  function  f(x,y )  is  to  be  determined  in  a  domain  [a  <  x  < 
6,  c  <  y  <  d]  from  the  knowledge  of  its  Fourier  coefficients  at  the  non-intersecting 
sub-domains  [at-  <  x  <  bi,Ci  <  y  <  d,}.  The  objective  is  to  approximate  the  “spliced” 
function  in  each  sub-domain  and  then  to  “glue”  the  approximations  together  in  order 
to  recover  the  original  function  in  the  full  domain. 

The  Fourier  partial  sum  approximation  in  each  sub-domain  yields  poor  results,  due  to 
the  Gibbs  phenomenon,  as  the  convergence  is  slow  and  spurious  oscillations  occur  at 
the  boundaries  of  each  sub-domain.  Thus  once  we  “glue”  the  sub-domain  approxima¬ 
tions  back  together,  the  approximation  for  the  function  in  the  full  domain  will  exhibit 
oscillations  throughout  the  entire  domain. 
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We  addressed  this  problem  in  one  and  two  dimensions  by  using  techniques  developed 
by  us  to  resolve  the  Gibbs  phenomenon.  We  have  shown  that  we  can  ’’cure”  the 
problem  completely  for  experimental  data.  The  results  of  this  effort  are  summarized 
in  [18] 

•  Parallel  Computing 

A  fast  direct  solver  has  been  developed  for  parallel  solution  of  “coarse  grid”  'problems, 
Aux  =  u/,,  such  as  arise  when  domain  decomposition  or  multi-grid  methods  are  ap¬ 
plied  to  elliptic  partial  differential  equations  in  d  space  dimensions.  The  approach  is 
based  upon  a  (quasi-)  sparse  factorization  of  the  inverse  of  A.  If  A  is  n  x  n  and  the 
number  of  processors  is  P,  our  approach  requires  0(n7  log2  P)  time  for  communica¬ 
tion  and  0(n1+7/P)  time  for  computation,  where  7  =  Results  from  a  512  node 
Intel  Paragon  show  that  our  algorithm  compares  favorably  to  more  commonly  used  ap¬ 
proaches  which  require  0(n  log2  P)  time  for  communication  and  0(n1+7)  or  0(n2/P) 
time  for  computation.  Moreover,  for  leading  edge  multicomputer  systems  with  thou¬ 
sands  of  processors  and  n  =  P  (i.e.,  communication  dominated  solves),  we  expect  our 
algorithm  to  be  markedly  superior  as  it  achieves  substantially  reduced  message  volume 
and  arithmetic  complexity  over  competing  methods  while  retaining  minimal  message 
startup  cost.  This  research  is  reported  in  [45]. 

Efficient  solution  of  the  Navier-Stokes  equations  in  complex  domains  is  dependent  upon 
the  availability  of  fast  solvers  for  sparse  linear  systems.  For  unsteady  incompressible 
flows,  the  pressure  operator  is  the  leading  contributor  to  stiffness,  as  the  characteristic 
propagation  speed  is  infinite.  In  the  context  of  operator  splitting  formulations,  it  is 
the  pressure  solve  which  is  the  most  computationally  challenging,  despite  its  elliptic 
origins. 

In  [14]  several  preconditioners  for  the  consistent  L2  Poisson  operator  arising  in  the 
spectral  element  formulation  of  the  incompressible  Navier-Stokes  equations  were  ex¬ 
amined.  A  finite  element  based  additive  Schwarz  preconditioner  using  overlapping 
sub-domains  plus  a  coarse  grid  projection  operator  which  is  applied  directly  to  the 
pressure  on  the  interior  Gauss  points,  was  developed.  For  large  two-dimensional  prob¬ 
lems  this  approach  can  yield  as  much  as  a  five-fold  reduction  in  simulation  time  over 
previously  employed  methods  based  upon  deflation. 

As  the  sound  speed  is  infinite  for  incompressible  flows,  computation  of  the  pressure 
constitutes  the  stiffest  component  in  the  time  advancement  of  unsteady  simulations. 
For  complex  geometries,  efficient  solution  is  dependent  upon  the  availability  of  fast 
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solvers  for  sparse  linear  systems.  In  [16]  a  Schwarz  preconditioner  for  the  spectral 
element  method  is  developed,  using  overlapping  sub-domains  for  the  pressure.  These 
local  sub-domain  problems  are  derived  from  tensor  products  of  one-dimensional  finite 
element  discretizations  and  admit  use  of  fast  diagonalization  methods  based  upon 
matrix-matrix  products.  In  addition,  we  use  a  coarse  grid  projection  operator  whose 
solution  is  computed  via  a  fast  parallel  direct  solver.  The  combination  of  overlapping 
Schwarz  preconditioning  and  fast  coarse  grid  solver  provides  as  much  as’  a  fourfold 
reduction  in  simulation  time  over  previously  employed  methods  based  upon  deflation 
for  parallel  architectures. 

•  General  Research 

It  is  the  nature  of  research  that  results  in  one  field  can  be  obtained  by  research  in 
different  fields.  Our  research  led  to  results  in  two  areas  not  intended  in  the  proposal. 

The  Nonlinear  Galerkin  Method  (NGM)  was  proposed  by  Temam  and  was  proven  to 
be  a  powerful  tool  in  approximating  complicated  dissipative  evolution  equations. 

The  method  is  based  on  the  idea  that  these  equations  describe  the  motion  of  small  scale 
as  well  as  those  of  large  scales.  The  method  suggest  the  factorization  of  the  equations 
to  those  describing  the  small  scales  and  those  of  the  large  scale.  The  equation  for  the 
small  scale  are  treated  differently,  since  not  the  small  scales  themselves  are  important 
but  rather  their  influence  on  the  large  scales. 

It  is  important  for  the  NGM  that  it  can  be  cast  within  the  framework  of  Spectral 
Methods.  This  had  been  done  in  joint  works  of  Gottlieb  and  Temam. 

In  [8]  it  had  been  shown  that  the  NGM  factorization  can  be  made  extremely  efficient 
if  different  time  advancing  methods  are  used  for  the  small  scales.  Several  explicit 
methods  are  suggested  for  the  small  scales  and  the  savings  are  outlined. 

In  [21]  the  super-convergence  of  the  Galerkin  methods  for  hyperbolic  initial  boundary 
value  PDE’s  has  been  discussed.  It  was  shown  that  super- convergence  is  lost  as  a 
result  of  the  imposition  of  boundary  conditions.  It  was  also  shown  that  there  is  no 
way  to  recover  the  super-convergence!!! 

In  [19]  the  shape  of  the  sea  surface  in  the  steady  state  solution  for  a  long  and  narrow 

•  basin,  as  the  Gulf  of  Suez  or  Baja  California,  is  studied.  The  study  addressed  the  time- 
dependent  problem  encountered  when  the  wind  in  the  wind  set-down  suddenly  relaxes 
and  the  water  gushes  landward  under  the  influence  of  the  pressure  gradient  force.  This 
is  a  difficult  problem  due  to  a  moving  singularity  associated  with  the  location  of  the 
intersection  point  between  the  sea  surface  and  the  shaping  bottom.  Spectral  methods 
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as  well  as  finite  difference  schemes  were  used.  Both  codes  yielded  the  same  results,  the 
spectral  methods  with  10  points  yielding  better  results  that  the  MacCormack  scheme 
with  3200  points!!  The  results  indicate  that  no  wave  breaking  occurs. 
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5.  Technology  Transfers 

•  Gottlieb  and  Don  cooperate  closely  with  A.  Nejad  and  Jeff  White  from  WPAFB  on 
the  problem  of  enhancing  mixing  by  shockwaves.  The  work  continues  with  Dr.  T. 
jackson. 

•  A  two  dimensional  code  for  subsonic  reactive  mixing  layer  was  delivered  for  Nejad’s 
group.  The  code  was  written  for  both  Cray  and  IBM  SP2  using  MPI.  This  was  a 
cooperative  effort  between  Nejad,  Givi  (Buffalo)  and  our  group. 

•  The  code  SPARC3D  is  used  in  Wright-Patterson  Lab.  It  is  an  extension  of  our  work 
on  fourth  order  schemes. 

9  We  have  started  a  joint  project  with  T.  Jackson  from  WPAFB  to  study  recessed  cavity 
fiameholders. 

•  We  have  close  cooperation  with  J.  Shang  concerning  the  construction  and  application 
of  compact  high  order  schemes. 
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•  We  have  started  a  cooperation  with  Young-Nam  Kim  from  Pratt  and  Whitney,  to 
implement  our  PML  method  for  acoustics  in  modeling  of  turbine  flutter. 
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